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ABSTRACT 

Latent  heat  energy  storage  systems  have  higher  energy 
density  than  their  sensible  heat  counterparts  and  have  the  added 
benefit  of  constant  temperature  operation.  This  work 
computationally  evaluates  a  thermal  energy  storage  system 
using  molten  silicon  as  a  phase  change  material.  A  cylindrical 
receiver,  absorber,  converter  system  was  evaluated  using  the 
heat  transfer  in  solids  with  surface-to  surface  radiation  physics 
module  of  the  commercially  available  COMSOL  Multiphysics 
simulation  software.  The  progression  of  the  solidification  and 
melting  fronts  through  the  phase  change  material  was  modeled 
for  two  different  methods  of  concentrated  solar  irradiation 
delivery.  Heating  the  core  of  the  PCM  rather  than  the  top  of  the 
PCM  decreased  the  required  solar  input  by  17%,  decreasing  the 
solar  collector  area  required  as  well  as  lowering  overall  system 
weight. 

INTRODUCTION 

In  order  to  use  solar  energy  as  a  constant  source  of  power, 
an  energy  storage  system  must  be  employed.  Thermal  energy 
storage  (TES)  allows  heat  energy  to  be  gathered  and  stored 
when  available,  and  utilized,  directly  or  converted  to  another 
form,  when  the  source  is  no  longer  available.  Historically, 
thermal  energy  has  been  stored  as  sensible  heat  in  either  water 
or  economical  solids,  such  as  rock  or  concrete  (1,  2).  Hasnain 
(2)  explains  that  a  sensible  heat  system  must  either  increase  the 
range  of  operating  temperatures  or  increase  the  mass  to  store 
more  energy,  thereby  limiting  its  use  in  a  small  system.  By 
comparison,  phase  change  materials  (PCMs)  utilize  the 
significant  heat  of  fusion  required  to  change  the  state  of  a 
material  (typically  solid  to  liquid);  this  allows  higher  energy 
storage  density  as  well  as  energy  storage  within  a  smaller 
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temperature  range.  In  2005,  the  International  Energy  Agency 
(IEA)  published  an  inventory  of  phase  change  materials 
appropriate  for  various  applications  including  energy  storage 

(3). 

The  use  of  PCM  in  low  temperature  energy  storage  for 
promoting  cost  effectiveness  and  energy  efficiency  in  buildings 
has  been  evaluated  (1,  2,  4).  However,  PCM  use  has  had  limited 
application  to  such  power  systems  due  to  materials  with  an 
appropriate  phase-change  temperature  (typically  hydrated  salts) 
offering  a  low  power  density  and  a  low  thermal  conductivity, 
leading  to  a  limited  rate  of  charging  and  discharging  (4).  A 
focus  on  developing  techniques  to  improve  PCM  conductivity 
has  resulted.  Imbedding  honeycomb  structures  in  the  PCM, 
macro-encapsulation,  and  microencapsulation  with  conductive 
materials  (1,  2,  3)  have  all  been  evaluated,  but  currently  the 
technology  requires  large  upfront  costs  for  small  increases  in 
efficiency  (4). 

One  field  with  perhaps  the  greatest  potential  for  successful 
application  of  PCM  technology,  however,  is  concentrated  solar 
power.  Pistocchini  (5)  investigated  the  use  of  PCM  storage  in 
the  field  of  concentrated  solar  power  plants.  In  this  terrestrial- 
based  application,  energy  can  be  stored  during  the  day  and 
exhausted  at  night,  capitalizing  on  a  large  desert  diurnal 
temperature  difference.  Similarly,  extraterrestrial  concentrated 
solar  power  is  another  application  area  with  the  potential  for 
significant  success.  For  small  satellites  in  low  earth  orbit,  where 
solar  exposure  is  followed  by  a  period  in  eclipse,  the  energy 
density  of  a  storage  system  is  paramount  for  launch  weight  (and 
therefore  cost)  reduction.  The  capabilities  of  solar  thermal 
energy  storage  in  the  form  of  a  PCM  with  high  melting 
temperatures  and  large  heats  of  fusion,  such  as  silicon,  could 
meet  these  needs  (6).  The  heat  of  fusion  and  thermal 
1  Copyright  ©  20 1 2  by  ASME 
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conductivity  of  silicon  are  one  and  two  orders  of  magnitude, 
respectively,  greater  than  those  of  hydrated  salts  and  paraffin 
waxes.  In  addition,  silicon  is  an  element  and  therefore  may 
avoid  the  material  degradation  experienced  with  other  phase 
change  materials.  The  material  properties  of  silicon  indicate 
great  potential  as  a  PCM,  warranting  a  more  detailed  numerical 
investigation. 

In  numerically  modeling  a  phase-changing  system, 
calculating  the  rate  of  movement  of  the  boundary  between  the 
melted  and  solid  portion  of  a  material  undergoing  phase  change 
is  referred  to  as  the  Stefan  problem.  Many  numerical 
techniques  have  been  applied  to  various  applications  of  the 
Stefan  problem,  including  the  relatively  simple  and  accurate 
heat  balance  integral  (HBIM)  (7)  and  enthalpy  methods  (8). 
These  two  techniques  were  applied  to  the  inward  and  outward 
cylindrical  solidification  case  by  Caldwell  and  Chan  (9)  and 
shown  to  be  consistent.  Ogoh  and  Groulx  presented  a  one¬ 
dimensional  (10)  and  a  cylindrical  finned  (11)  model  of  the 
phase  change  process  in  the  commercially  available  software 
COMSOL.  To  simulate  phase  change  they  used  the  heat 
capacity  method.  Their  one-dimensional  results  compared  well 
to  an  analytical  evaluation. 

The  computational  work  described  here  simulates  a  PCM- 
based  thermal  storage  system  using  silicon  in  low  earth  orbit. 
An  altitude  of  1900  km  was  selected  to  represent  the  most 
demanding  case  of  the  low  earth  orbit  with  the  longest  time 
spent  in  eclipse.  At  this  altitude  4500  seconds  are  spent  in 
daylight  and  3000  seconds  in  eclipse.  The  progression  of  the 
solidification  and  melting  fronts  through  the  phase  change 
material  was  modeled  to  determine  the  amount  of  latent  heat 
energy  remaining  in  the  core  and  subsequently  the  required 
energy  input. 

MODEL 

The  design  of  a  system  capable  of  receiving,  absorbing, 
and  converting  solar  energy,  referred  to  as  the  RAC,  was 
evaluated  computationally.  Describing  the  system  from  the 
center  out  radially,  a  47  mm  radius  rod  of  phase  change 
material  with  a  height  of  92  mm  was  housed  in  a  cylindrical 
inner  container  with  walls  5  mm  thick.  A  25  mm  vacuum  gap 
separated  the  outer  walls  of  the  inner  container  from  two 
radiation  shields,  a  55  mm  thick  layer  of  insulation,  and  an 
outer  container  with  10  mm  thick  walls.  A  thermal  to  electric 
conversion  (TEC)  system  were  located  above  the  inner 
container  top  surface.  The  TEC  device  would  be  a  low 
bandwidth  photovoltaic  cell  with  a  band  gap  of  approximately 
0.55eV,  best  simulated  by  an  InGaAsSb  cell.  A  top-down  view 
of  the  system  is  shown  in  Figure  1  and  a  side  view  in  Figure  3. 
The  phase  change  material  was  silicon,  the  inner  container 
silicon  carbide,  the  radiation  shields  gold,  the  insulation 
carbon-bonded  carbon  fiber,  and  the  outer  casing  graphite.  The 
system  was  designed  so  radial  energy  losses  were  minimized. 

A  two-dimensional  model,  generated  in  the  commercially 
available  COMSOL  Multiphysics,  was  used  to  evaluate  the 
melting  and  solidification  fronts  in  the  phase  change  material 
and  irradiation  of  the  thermal  to  electric  conversion  system. 


One-half  of  the  test  section  cross-section  was  modeled  and  an 
axial  symmetry  condition  applied  to  the  centerline.  A  surface- 
to-surface  radiation  condition  was  applied  to  the  outer 
circumference  of  the  inner  container  of  the  system,  all  surfaces 
of  both  radiation  shields,  the  bottom  surface  of  the  TEC,  and 
the  inner  circumference  of  the  insulation  layer.  A  surface  to 
ambient  radiation  condition  was  applied  to  the  outer 
circumference  of  the  outer  container  as  well  as  its  top  and 
bottom  and  the  top  of  the  TEC.  Incoming  concentrated  solar 
energy  was  modeled  using  two  different  methods,  a  constant 
flux  was  applied  to  the  top  of  the  inner  container  for  Method  I 
and  a  constant  flux  applied  on  the  bottom  of  a  narrow 
cylindrical  extension  of  the  container  which  extends  into  the 
center  of  the  PCM  for  Method  II.  For  all  conditions,  an  ambient 
temperature  of  300  K  was  assigned. 

A  combination  of  three  different  studies,  one  steady  and 
two  transient,  all  using  the  COMSOL  heat  transfer  in  solids 
physics  module,  were  used  to  model  the  phase  change  process. 
A  steady  state  model  where  a  constant  temperature  condition 
was  applied  to  the  outer  surface  of  the  phase  change  material 
was  evaluated  to  determine  temperature  profiles  in  the  system 
after  melting.  These  temperature  profiles  were  used  as  an  input 
to  the  transient  study.  In  the  transient  study  the  constant  surface 
temperature  setting  of  the  steady  state  evaluation  was  removed 
and  the  system  was  allowed  to  cool  for  the  3000  s  associated 
with  eclipse.  The  temperature  profile  after  the  completion  of 
this  phase  was  used  as  the  initial  condition  for  the  heating 
period.  The  heat  was  applied  as  a  constant  heat  flux,  in  one  of 
the  two  methods  described  above,  for  the  duration  of  the  4500  s 
daylight  phase. 

To  mitigate  the  effects  of  beginning  the  orbit  in  the  ideal 
case  of  the  steady  state  solution  of  a  perfectly  melted  core,  a 
second  orbit  was  simulated  with  the  final  profile  from  the 
heating  cycle  serving  as  the  initial  condition  for  the  cooling 
phase  and  that  fully  cooled  profile  subsequently  being  used  as 
the  initial  condition  for  a  second  heating  phase.  A  direct 
PARDISO  solver  was  used  for  all  the  components  of  the  study. 

The  heat  capacity  method,  where  the  specific  heat  of  the 
liquid  material,  CPji,  is  increased  by  the  latent  heat  of  fusion,  lf, 
over  a  temperature  range,  A Tm,  centered  on  the  phase  change 
temperature  was  used  to  simulate  phase  change. 


LA1m 

When  the  temperature  of  a  particular  element  was  within  half 
A Tm  of  the  melting  temperature  a  step  discontinuity  in  the 
specific  heat  occurred.  The  increase  in  specific  heat  was 
modeled  using  a  Heaviside  smoothed  step  function  over  a  1 0  K 
temperature  range. 

VALIDATION 

The  steady  state  heat  flux  through  different  insulation  and 
insulation/radiation  shield  configurations  was  evaluated  in 
MATLAB  and  COMSOL  for  an  initial  validation  of  the 
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radiation  modeling  in  COMSOL  as  well  as  to  determine  a 
design  configuration  that  minimized  thermal  losses  to  the 
environment.  The  design  described  in  the  model  section  is  the 
result  of  this  preliminary  validation. 

The  system  was  modeled  as  concentric,  one-dimensional 
cylinders.  In  all  cases,  the  inner  cylinder  surface  temperature 
was  set  to  1000  K  and  the  outer  surface  radiated  to  an  ambient 
temperature  of  0  K.  Emissivities  and  thermal  conductivities 
were  selected  to  model  the  design  materials  and  available 
coatings.  The  outer  surface  temperature  and  heat  flux  were 
predicted.  Several  different  configurations  were  evaluated. 

The  best  results  were  obtained  by  insulating  the  hot  side  of 
the  vacuum  gap  and  utilizing  polishing  or  coating  to  reduce 
radiative  heat  transfer,  however,  the  high  temperature  of  the 
container  in  this  case  made  this  configuration  unusable.  Not 
insulating  the  outer  circumference  of  the  inner  container,  and 
instead  using  radiation  shielding  between  the  container  and  the 
insulation,  resulted  in  only  slightly  larger  thermal  losses.  The 
results  of  this  brief  study  lead  to  the  combination  of  insulation 
on  the  outer  casing  to  avoid  exposure  to  peak  temperatures  and 
multiple  radiation  shields.  Comparison  of  the  analytical  and 
computational  predictions  show  good  agreement  with 
differences  in  temperature  predictions  less  than  0.35%  and  in 
the  predicted  energy  loss  of  1.32%. 

Inner  container  Radiation 


Insulation  Vacuum  gap  Insulation 


FIGURE  1:  Top  down  view  of  (a)  best  performing  and  (b) 
final  design  configurations. 

Once  the  steady  state  radiation  model  was  validated,  the 
phase  change  model  employed  in  COMSOL  was  evaluated.  The 
solidification  of  a  cylinder  of  radius,  r0,  exposed  to  a  constant 
surface  temperature  was  simulated  in  COMSOL.  The  location 
of  the  solidification  front  as  a  function  of  time  was  predicted 
and  compared  to  the  enthalpy  model  for  inward  cylindrical 
solidification  published  by  Caldwell  and  Chen  (9)  which 
showed  excellent  agreement  with  the  cylindrical  HBIM  model. 
The  COMSOL  results  are  compared  with  those  of  Caldwell  and 
Chen  in  Figure  2.  The  COMSOL  results  agree  well  with  the 
numerical  values,  validating  the  phase  change  portion  of  the 
model. 


0  0.1  0.2  0.3  0.4 

Non-Dimensional  Time 


FIGURE  2:  Comparison  of  solidifcation  front  predicted  by 
Caldwell  (9)  and  COMSOL  model. 

RESULTS  AND  DISCUSSION 

A  receiver,  absorber,  converter  system  in  low  earth  orbit 
was  modeled.  An  altitude  of  1900  km  was  selected,  resulting  in 
4500  s  spent  in  daylight  and  3000  s  in  eclipse  during  each  orbit. 
The  progression  of  the  solidification  and  melting  fronts  through 
the  phase  change  material  during  orbit  was  modeled  to  track  the 
amount  of  latent  heat  energy  remaining  in  the  core.  The  results 
from  the  cooling  and  heating  cycles  are  presented  for  two 
different  cases  (1)  top  heating  and  (2)  center  heating.  In  the  top 
heating  case,  Method  I,  a  concentrated  heat  flux  is  applied  to 
the  top  of  the  inner  container.  In  the  center  heating  case, 
Method  II,  a  heat  flux  is  applied  to  the  center  of  the  PCM. 

Method  I:  Top  Heating 

In  the  top  heating  study,  the  heat  flux  applied  to  the  inner 
container  top  was  varied  until  the  PCM  melted  completely 
during  the  daylight  phase  and  fully  solidified  during  the  eclipse 
phase.  This  corresponded  to  a  constant  surface  heat  transfer  rate 
of  1800  W.  The  results  presented  are  those  found  with  the  1800 
W  input. 

The  temperature  profile  along  the  cutline,  shown  in  Figure 
3,  associated  with  solidification  during  eclipse  is  shown  in 
Figure  4  as  a  function  of  time.  Energy  is  lost  from  the  inner 
container  circumference,  top,  and  bottom  through  radiation  and 
minor  conduction  through  the  standoffs  by  which  the  core  is 
suspended.  Because  the  system  is  heated  from  the  top,  as 
solidification  begins,  the  top  of  the  system  is  significantly 
warmer.  The  temperature  of  the  center  of  the  silicon  core 
reaches  the  phase  change  temperature  in  less  than  300  s  but 
requires  the  remainder  of  the  3000  s  phase  to  completely 
solidify. 
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FIGURE  3:  Side  view  of  RAC  system  with  cut  line  used  for 
data  collection  identified. 

The  progression  of  the  solidification  front  from  the  center  out  is 
shown  in  Figure  5.  The  different  shades  show  an  isotherm 
associated  with  1687  K  at  different  instances  of  time.  The 
solidification  front  progresses  more  quickly  from  the  top 
because  the  placement  of  the  TEC  prevents  the  installation  of 
insulation  and  radiation  shielding  at  the  top  of  the  container. 


Line  Graph:  Temper  alure  (K) 


FIGURE  4:  Temperature  profile  along  the  cut  line 
identified  in  Figure  3  at  different  times  during 
solidification. 
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FIGURE  5:  Progression  of  solidification  front  with  time. 

The  solidification  and  sensible  heat  loss  from  the  system 
result  in  the  temperature  profile  shown  in  Figure  4  associated 
with  time  equal  to  3000  s.  This  temperature  profile  serves  as  the 
initial  condition  for  the  daylight  phase  of  the  orbit.  The 
temperature  profile  along  the  cutline  through  the  PCM  only, 
shown  in  Figure  6,  shows  the  change  in  the  temperature  profile 
as  the  top  surface  of  the  inner  container  is  subjected  to  a 
constant  surface  heat  flux  during  the  4500  s  of  solar  irradiation. 

The  thermal  resistance  of  silicon  leads  to  a  large 
temperature  gradient  between  the  top  and  bottom  of  the  molten 
core.  The  temperature  at  the  interface  between  the  casing  and 
the  molten  core  is  approaching  the  melting  temperature  of  the 
casing,  which  threatens  the  structure  of  the  device. 


FIGURE  6:  Temperature  profile  along  the  cut  line  through 
the  PCM  identified  in  Figure  3  at  different  times  during 
melt. 
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The  melting  front  is  shown  with  isotherms  associated  with 
1687  K  in  Figure  7.  The  phase  change  front  moves  from  the 
inner  container  top  vertically  towards  the  bottom.  As  the 
melting  front  moves  down,  the  front  moves  more  quickly  along 
the  edges  than  along  the  centerline  because  the  silicon  carbide 
casing  is  approximately  six  times  more  conductive  than  the 
silicon  core. 


FIGURE  7:  Progression  of  solidification  front  with  time. 


The  heat  added  to  the  top  surface  of  the  container  both  conducts 
downward  into  the  PCM  but  also  through  the  standoff.  In 
addition,  the  standoff  is  subjected  to  radiation  from  the  top  of 
the  container.  The  combined  conduction  and  radiation  raises  the 
standoff  temperature  well  above  the  material  melting  point  as  is 
indicated  by  the  shading  in  the  standoffs  in  Figure  7.  Because 
this  analysis  showed  many  of  the  system  materials  reaching 
their  melting  points,  it  was  determined  that  a  different  heating 
method  would  be  required,  leading  to  the  center  heating  study 
described  below. 

Method  II:  Center  Heating 

In  the  center  heating  study,  the  heat  flux  applied  to  the 
center  of  the  system  was  varied  until  the  PCM  melted 
completely  during  the  daylight  phase  and  hilly  solidified  during 
the  eclipse  phase.  This  corresponded  to  a  constant  surface  heat 
transfer  rate  of  1500  W.  The  results  presented  are  those  found 
with  the  1500  W  input. 

The  center  heating  method  concentrates  the  solar  heat  flux 
on  the  base  of  a  vertical  channel  located  in  the  center  of  the 
PCM.  The  temperature  profile  along  the  cutline  in  the  PCM, 
shown  in  Figure  8,  associated  with  solidification  of  the  center 
heated  PCM  during  eclipse  is  shown  in  Figure  9  as  a  function 
of  time.  The  temperature  profile  at  time  equal  to  zero  seconds  is 
the  final  temperature  profile  associated  with  the  previous  melt 
phase.  After  approximately  300  s,  the  entire  cutline  in  the  core 
is  at  the  melt  temperature  and  at  each  time  step  afterward,  the 
solidification  front  moves  downward.  The  temperatures  are 


again  lower  on  the  top  of  the  core  than  the  bottom  as  the  TEC 
system  prevents  the  installation  of  insulation  and  allows  for 
losses. 


FIGURE  8:  Side  view  of  RAC  system  with  cut  line  used  for 
data  collection  identified. 


FIGURE  9:  Temperature  profile  along  the  cut  line 
identified  in  Figure  8  at  different  times  during 
solidification. 
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I  Ine  Graph:  Temperature  (K) 


FIGURE  10:  Temperature  profile  along  the  cut  line 
identified  in  Figure  8  at  different  times  during  melt. 


The  temperature  profiles  during  the  melt,  as  shown  in 
Figure  10,  increase  symmetrically  with  the  highest 
temperatures  at  the  PCM  center.  Melting  from  the  center  out  is 
beneficial  as  the  solid  silicon  acts  as  insulation  during  the 
melting  process  and  interfacial  issues  are  eliminated.  In 
addition,  the  required  heat  flux  is  reduced  by  about  17%  which 
correspondingly  decreases  the  required  collector  area  and 
associated  weight. 

Figure  11  shows  the  expanding  melting  front  during  the 
center  heating  method.  The  different  shades  show  a  1687  K 
isotherm  at  different  times.  Excessive  temperatures  do  not 
reach  the  standoff,  but  the  majority  of  the  PCM  completely 
melts  and  solidifies  taking  advantage  of  the  latent  heat  storage 
without  risking  melting  of  any  structural  elements. 

Melting  cycle 


FIGURE  11:  Progression  of  melt  front  with  time. 

The  results  for  the  center  heated  configuration  predicts  a 
container  top  temperture  that  oscillates  between  1550K  and 
1650K  during  a  single  orbit.  Therefore,  the  TEC  surface 


receives  an  average  2000  W  of  irradiation  for  this  case, 
resulting  in  an  approximate  electric  power  output  of  200  W 
based  on  a  50%  spectral  efficiency  and  a  20%  TPV  cell 
efficiency  resulting  in  a  conservative  estimate  of  10%  electric 
conversion  efficiency.  In  order  to  compare  this  system  to  a 
comparable  design  powered  by  a  photovoltaic  array  and  bank 
of  batteries,  a  brief  analysis  was  conducted  to  determine  the 
required  solar  panel  size  and  system  weight.  In  the  same  orbit,  a 
3  m2  array  of  photovoltaic  cells  would  be  required  that  in 
combination  with  battery  storage  would  weigh  approximately 
20  kg.  By  comparison,  the  silicon  design  with  a  70%  solar 
concentration  efficiency  (12)  would  require  2.2  m2  of  collector 
area  and  a  total  weight  of  approximately  1 8  kg. 

CONCLUSION 

A  computational  model  of  a  thermal  energy  storage  system 
using  molten  silicon  as  the  phase  change  material  was 
presented.  The  system  was  modeled  in  low  earth  orbit  with 
4500  s  spent  in  daylight  and  3000  s  in  eclipse.  The  progression 
of  the  solidification  and  melting  fronts  through  the  phase 
change  material  was  modeled  for  two  different  methods  of 
concentrated  solar  irradiation  delivery.  Heating  the  core  of  the 
PCM  rather  than  the  top  of  the  PCM  decreased  the  required 
solar  input  by  17%  which  would  decrease  the  solar  collector 
area  required  and  lower  overall  system  weight.  This 
preliminary  study  will  be  expanded  to  include  convection  in  the 
liquid  phase  as  well  as  to  investigate  more  suitable  insulation 
materials. 

NOMENCLATURE 

Cp  Specific  heat  (J/kg/K) 

L  Latent  heat  (J/kg) 

T  Temperature  (K) 

Superscripts  and  Subscripts 
f  Fusion 

1  Liquid 

m  Melt 

Greek  Letters 
A  Change 
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